package cn.az13js.satool;

import java.io.File;
import java.util.LinkedList;
import java.io.IOException;

public class App {

    public static void main(String[] args) throws IOException {
        if (Util.commandParamExists(args, "-test_debug")) {
            double[] a = Util.nDArray(100);
            double[] x = {0.5, -0.2, 0.3};
            for (int i = 0; i < a.length; i++) {
                a[i] = Util.normalDistribution(0, 0.1);
                for (int j = 0; j < x.length; j++) {
                    if (i - j - 1 > -1) {
                        a[i] += x[j] * a[i - j - 1];
                    }
                }
                Util.message("" + a[i]);
            }
            return;
        }
        if (0 == args.length || Util.commandParamExists(args, "-h")) {
            printHelpMessage();
            return;
        }
        if (!Util.commandParamExists(args, "-i") || !Util.commandParamExists(args, "-o")) {
            printHelpMessage();
            return;
        }
        LinkedList<Double> dataFromCsv = Util.readData(Util.commandParam(args, "-i"));
        double[] dataFromCsvArray = new double[dataFromCsv.size()];
        int i = 0;
        for (double a : dataFromCsv) {
            dataFromCsvArray[i++] = a;
        }
        int free = -1;
        if (Util.commandParamExists(args, "-m")) {
            free = Integer.parseInt(Util.commandParam(args, "-m"));
            if (free >= dataFromCsvArray.length || free < 1) {
                Util.message("“-m”指定的自由度范围超出限制，自由度只能取大于0和小于序列长度的值。");
                free = -1;
            }
        }
        int d = 200;
        if (Util.commandParamExists(args, "-d")) {
            d = Integer.parseInt(Util.commandParam(args, "-d"));
            if (d < 1) {
                Util.message("“-d”指定的概率分布精细度超出限制，其只能取大于0。");
                d = 200;
            }
        }
        run(dataFromCsvArray, Util.commandParam(args, "-o"), free, d);
    }

    public static double run(double[] seq, String outputDir, int free, int dis) throws IOException {
        String sequenceFileName = outputDir + File.separator + "sequence.csv"; // 存储序列本身的文件名
        String distributionFileName = outputDir + File.separator + "distribution.csv"; // 存储数据分布的文件名
        String afFileName = outputDir + File.separator + "af.csv"; // 存储数据AF的文件名
        String acfFileName = outputDir + File.separator + "acf.csv"; // 存储数据ACF的文件名
        String pacfFileName = outputDir + File.separator + "pacf.csv"; // 存储数据PACF的文件名
        String qFileName = outputDir + File.separator + "q.csv"; // 存储数据卡方分析的文件名
        // 测试序列长度
        int len = seq.length;
        double minValue = seq[0];
        double maxValue = seq[0];
        for (double a : seq) {
            if (a < minValue) {
                minValue = a;
            }
            if (a > maxValue) {
                maxValue = a;
            }
        }
        Util.message("\"********** Basic info **********\"");
        Util.message("LEN=" + len);
        Util.message("MIN=" + minValue);
        Util.message("MAX=" + maxValue);
        // 计算平均值
        double avg = 0;
        for (double e : seq) {
            avg += e;
        }
        avg /= len;
        Util.message("AVG=" + avg);
        // 计算方差
        double var = 0;
        for (double e : seq) {
            var += (e - avg) * (e - avg);
        }
        var /= len;
        Util.message("VAR=" + var);
        double[][] sequenceArray = Util.nDArray(len, 2);
        for (int i = 0; i < len; i++) {
            sequenceArray[i][0] = i + 1;
            sequenceArray[i][1] = seq[i];
        }
        WriteACSVFile.write(sequenceFileName, sequenceArray);
        WriteACSVFile.write(distributionFileName, Distribution.getAnalyzeInfo(seq, minValue, maxValue, dis));
        // 计算AF和ACF
        double[] af = Util.nDArray(len);
        double[][] af2std = Util.nDArray(len, 2); // 0是负2倍标准差偏移，1是正2倍标准差偏移
        double[] acf = Util.nDArray(len);
        double[][] acf2std = Util.nDArray(len, 2); // 0是负2倍标准差偏移，1是正2倍标准差偏移
        af[0] = var;
        af2std[0][0] = var - Math.sqrt(2 * var * var / len);
        af2std[0][1] = var + Math.sqrt(2 * var * var / len);
        acf[0] = 1;
        acf2std[0][0] = 1;
        acf2std[0][1] = 1;
        double avg1 = 0;
        double avg2 = 0;
        double var1 = 0;
        double var2 = 0;
        double sum = 0;
        double temp = 0;
        for (int k = 1; k < len; k++) {
            avg1 = 0;
            avg2 = 0;
            sum = 0;
            for (int i = 0; i < len - k; i++) {
                avg1 += seq[i];
                avg2 += seq[i + k];
                sum += seq[i] * seq[i + k];
            }
            avg1 = avg1 / (len - k);
            avg2 = avg2 / (len - k);
            sum = sum / (len - k);
            af[k] = sum - avg1 * avg2;
            if (k >= len - 1) {
                af2std[k][0] = 0;
                af2std[k][1] = 0;
            } else if (k > len / 2) {
                af2std[k][0] = -1 * Math.sqrt(var * var / 4);
                af2std[k][1] = Math.sqrt(var * var / 4);
            } else {
                af2std[k][0] = -1 * Math.sqrt(var * var / (len - k));
                af2std[k][1] = Math.sqrt(var * var / (len - k));
            }
            temp = Math.sqrt(1.0 / (double)(len - k + 2)); // 标准差
            if (k < len - 1) {
                var1 = 0;
                var2 = 0;
                for (int i = 0; i < len - k; i++) {
                    var1 += (seq[i] - avg1) * (seq[i] - avg1);
                    var2 += (seq[i + k] - avg2) * (seq[i + k] - avg2);
                }
                var1 /= len - k;
                var2 /= len - k;
                acf[k] = af[k] / Math.sqrt(var1 * var2);
                if (k <= len / 4) {
                    acf2std[k][0] = -1.0 / (double)len - 2.0 * temp;
                    acf2std[k][1] = -1.0 / (double)len + 2.0 * temp;
                } else if (k <= len / 2) {
                    acf2std[k][0] = (double)(4 * k - 2 * len) / len / len - 2.0 * temp;
                    acf2std[k][1] = (double)(4 * k - 2 * len) / len / len + 2.0 * temp;
                } else {
                    acf2std[k][0] = -2.0 * temp;
                    acf2std[k][1] = 2.0 * temp;
                }
            } else {
                acf[k] = 0;
                acf2std[k][0] = -2.0 * temp;
                acf2std[k][1] = 2.0 * temp;
            }
        }
        Util.message("\"********** Advanced **********\"");
        double[][] seqAf = Util.nDArray(af.length, 4);
        for (int i = 0; i < af.length; i++) {
            seqAf[i][0] = i;
            seqAf[i][1] = af[i];
            seqAf[i][2] = af2std[i][0];
            seqAf[i][3] = af2std[i][1];
        }
        WriteACSVFile.write(afFileName, seqAf);
        double[][] seqAcf = Util.nDArray(acf.length, 4);
        for (int i = 0; i < acf.length; i++) {
            seqAcf[i][0] = i;
            seqAcf[i][1] = acf[i];
            seqAcf[i][2] = acf2std[i][0];
            seqAcf[i][3] = acf2std[i][1];
        }
        WriteACSVFile.write(acfFileName, seqAcf);
        // 计算PACF
        double[][] pacf = new double[len - 1][];
        double[] error = new double[len - 1];
        for (int k = 0; k < len - 1; k++) {
            pacf[k] = new double[k + 1];
            if (0 == k) { // PACF(1), PACF(1) -> 程序double[][]的pacf[0][0]
                pacf[k][k] = af[1] / af[0];
                error[k] = af[0] * (1 - pacf[k][k] * pacf[k][k]);
            } else {
                temp = 0;
                for (int i = 0; i < k; i++) {
                    temp += af[k - i] * pacf[k - 1][i];
                }
                pacf[k][k] = (af[k + 1] - temp) / error[k - 1];
                for (int i = 0; i < k; i++) {
                    pacf[k][i] = pacf[k - 1][i] - pacf[k][k] * pacf[k - 1][k - i - 1];
                }
                error[k] = error[k - 1] * (1 - pacf[k][k] * pacf[k][k]);
            }
        }
        double[][] seqPacf = Util.nDArray(pacf.length, 5);
        for (int i = 0; i < pacf.length; i++) {
            seqPacf[i][0] = 1 + i;
            seqPacf[i][1] = pacf[i][i];
            seqPacf[i][2] = error[i] > 0 ? -2.0 * Math.sqrt(error[i]) : 0;
            seqPacf[i][3] = error[i] > 0 ? 2.0 * Math.sqrt(error[i]) : 0;
            seqPacf[i][4] = error[i];
        }
        WriteACSVFile.write(pacfFileName, seqPacf);
        // 自由度计算输出
        int m;
        if (free < 1) {
            m = (int)Math.ceil(Math.log(len));
        } else {
            m = free;
        }
        Util.message("m=" + m);
        temp = 0;
        for (int i = 1; i <= m; i++) {
            temp += (acf[i] * acf[i] / (double)(len - i));
        }
        double Q = len * (len + 2) * temp;
        double[][] seqQ = Util.nDArray(20, 3);
        seqQ[0][2] = 3.841;
        seqQ[1][2] = 5.991;
        seqQ[2][2] = 7.815;
        seqQ[3][2] = 9.488;
        seqQ[4][2] = 11.071;
        seqQ[5][2] = 12.592;
        seqQ[6][2] = 14.067;
        seqQ[7][2] = 15.507;
        seqQ[8][2] = 16.919;
        seqQ[9][2] = 18.307;
        seqQ[10][2] = 19.675;
        seqQ[11][2] = 21.026;
        seqQ[12][2] = 22.362;
        seqQ[13][2] = 23.685;
        seqQ[14][2] = 24.966;
        seqQ[15][2] = 26.296;
        seqQ[16][2] = 27.587;
        seqQ[17][2] = 28.869;
        seqQ[18][2] = 30.144;
        seqQ[19][2] = 31.410;
        double q;
        // 自动生成按m为自变量，Q为因变量的函数
        for (int i = 1; i <= 20; i++) {
            temp = 0;
            for (int j = 1; j <= i; j++) {
                temp += (acf[j] * acf[j] / (double)(len - j));
            }
            q = len * (len + 2) * temp;
            seqQ[i - 1][0] = i;
            seqQ[i - 1][1] = q;
        }
        WriteACSVFile.write(qFileName, seqQ);
        Util.message("Q=" + Q);
        return 0;
    }

    private static void printHelpMessage() {
        Util.message("命令行参数：\n-i 数据文件名（必填）\n-o 输出目录（必填）\n-h 打印此帮助信息\n-m 自由度，默认为长度的自然对数值向上取整\n-d 统计概率分布时横坐标的精度\n示例：program -i input.csv -o output_dir\n");
    }
}
